{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example: combining MOD file ion channels with rxd"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A version of this notebook may be run online via Google Colab at https://tinyurl.com/neuron-rxd-and-mod\n",
" (make a copy or open in playground mode)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Overview"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"NEURON's reaction-diffusion infrastructure can be used to readily allow intracellular concentrations to respond to currents generated in MOD files, as long as:\n",
"\n",
"- `nrn_region='i'` is specified for the `rxd.Region` (so that it knows it corresponds to the electrophysiology region of the inside of the cell), AND\n",
"the name and charge of the ion/etc are given in the `rxd.Species` declaration.\n",
"Satisfying the above two rules also allows MOD files to see intracellular concentrations.\n",
"\n",
"- 3D extracellular concentrations also interoperate with electrophysiology automatically as long as name and charge are specified.\n",
"\n",
"As a simple example, we consider a model with just a single point soma, of length and diameter 10 microns, with Hodgkin-Huxley kinetics (which uses the built in mod file `hh.mod`), and dynamic sodium (declared using rxd but without any additional kinetics)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup NEURON library and imports"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's import our usual NEURON libraries and definitions. Remember you can use either `um` or `µm` for micron."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2025-05-23T00:19:32.561974Z",
"iopub.status.busy": "2025-05-23T00:19:32.561830Z",
"iopub.status.idle": "2025-05-23T00:19:32.869507Z",
"shell.execute_reply": "2025-05-23T00:19:32.869144Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"1.0"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from neuron import h, rxd\n",
"from neuron.units import mV, ms, um, mM\n",
"\n",
"## needed for standard run system\n",
"h.load_file(\"stdrun.hoc\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now import `plotly`, a graphics library. (You could easily modify this code to use other graphics libraries like `matplotlib`, `plotnine`, or `bokeh`.)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2025-05-23T00:19:32.871231Z",
"iopub.status.busy": "2025-05-23T00:19:32.871031Z",
"iopub.status.idle": "2025-05-23T00:19:32.909417Z",
"shell.execute_reply": "2025-05-23T00:19:32.909052Z"
}
},
"outputs": [],
"source": [
"import plotly.graph_objects as go\n",
"from plotly.subplots import make_subplots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup the model"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2025-05-23T00:19:32.911322Z",
"iopub.status.busy": "2025-05-23T00:19:32.910934Z",
"iopub.status.idle": "2025-05-23T00:19:32.914908Z",
"shell.execute_reply": "2025-05-23T00:19:32.914588Z"
}
},
"outputs": [],
"source": [
"## define morphology\n",
"soma = h.Section(name=\"soma\")\n",
"soma.L = soma.diam = 10 * um\n",
"\n",
"## add ion channels (h.hh is built in, so always available)\n",
"h.hh.insert(soma)\n",
"\n",
"## define cytosol. MUST specify nrn_region for concentrations to update\n",
"cyt = rxd.Region([soma], name=\"cyt\", nrn_region=\"i\")\n",
"\n",
"## define sodium. MUST specify name and charge for concentrations to update\n",
"na = rxd.Species(cyt, name=\"na\", charge=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Alternatively, we could have written `h.hh.insert(soma.wholetree())` to put Hodgkin-Huxley channels everywhere in the cell that the soma is part of, but since there is only one section, this is not required."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's also add an excitatory synapse to receive events (these will trigger the cell to spike)."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2025-05-23T00:19:32.916339Z",
"iopub.status.busy": "2025-05-23T00:19:32.916164Z",
"iopub.status.idle": "2025-05-23T00:19:32.918446Z",
"shell.execute_reply": "2025-05-23T00:19:32.918071Z"
}
},
"outputs": [],
"source": [
"syn = h.ExpSyn(soma(0.5))\n",
"syn.tau = 1 * ms\n",
"syn.e = 0 * mV"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Add a stimulus"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The spike events themselves (two events, 15 ms apart starting at h.t=10*ms):"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2025-05-23T00:19:32.919829Z",
"iopub.status.busy": "2025-05-23T00:19:32.919543Z",
"iopub.status.idle": "2025-05-23T00:19:32.922114Z",
"shell.execute_reply": "2025-05-23T00:19:32.921798Z"
}
},
"outputs": [],
"source": [
"stim = h.NetStim()\n",
"stim.interval = 15 * ms\n",
"stim.number = 2\n",
"stim.start = 10 * ms"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Send those events to our synapse:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2025-05-23T00:19:32.923601Z",
"iopub.status.busy": "2025-05-23T00:19:32.923259Z",
"iopub.status.idle": "2025-05-23T00:19:32.926113Z",
"shell.execute_reply": "2025-05-23T00:19:32.925735Z"
}
},
"outputs": [],
"source": [
"nc = h.NetCon(stim, syn)\n",
"nc.weight[0] = 0.001"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup recording variables"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2025-05-23T00:19:32.928155Z",
"iopub.status.busy": "2025-05-23T00:19:32.927505Z",
"iopub.status.idle": "2025-05-23T00:19:32.930701Z",
"shell.execute_reply": "2025-05-23T00:19:32.930361Z"
}
},
"outputs": [],
"source": [
"t = h.Vector().record(h._ref_t)\n",
"v = h.Vector().record(soma(0.5)._ref_v)\n",
"na_vec = h.Vector().record(soma(0.5)._ref_nai)\n",
"ina = h.Vector().record(soma(0.5)._ref_ina)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Initialize and run the simulation"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"execution": {
"iopub.execute_input": "2025-05-23T00:19:32.932629Z",
"iopub.status.busy": "2025-05-23T00:19:32.931992Z",
"iopub.status.idle": "2025-05-23T00:19:32.940869Z",
"shell.execute_reply": "2025-05-23T00:19:32.940529Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"h.finitialize(-65 * mV)\n",
"h.continuerun(50 * ms)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plot it"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"execution": {
"iopub.execute_input": "2025-05-23T00:19:32.942946Z",
"iopub.status.busy": "2025-05-23T00:19:32.942311Z",
"iopub.status.idle": "2025-05-23T00:19:33.122485Z",
"shell.execute_reply": "2025-05-23T00:19:33.122159Z"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" \n",
" "
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig = make_subplots(rows=3, cols=1)\n",
"fig.add_trace(go.Scatter(x=t, y=v, name=\"v\"), row=1, col=1)\n",
"fig.update_yaxes(title_text=\"v (mV)\", row=1, col=1)\n",
"fig.add_trace(go.Scatter(x=t, y=ina, name=\"ina\"), row=2, col=1)\n",
"fig.update_yaxes(title_text=\"ina (mA/cm^2)\", row=2, col=1)\n",
"fig.add_trace(go.Scatter(x=t, y=na_vec, name=\"[Na+]\"), row=3, col=1)\n",
"fig.update_xaxes(title_text=\"t (ms)\", row=3, col=1)\n",
"fig.update_yaxes(title_text=\"[Na+] (mM)\", row=3, col=1)\n",
"\n",
"fig.show(renderer=\"notebook_connected\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Note: On homeostatic mechanisms or the lack thereof\n",
"\n",
"Without any additional homeostatic mechanisms (Hodgkin and Huxley did not model sodium concentration so they did not need to include homeostatic mechanisms for it), intracellular sodium concentration will not return to baseline, and each spike will move intracellular sodium concentration closer to the extracellular concentration. Potassium concentration in this model is constant as we did not declare a potassium rxd.Species, but if we did it would also approach its extracellular concentration with each spike, until eventually the cell is not able to fire action potentials anymore."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Finally: The mod file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For any distributed mechanism, including the built-in `h.hh`, one can always get the source code for the corresponding mod file using the `.code` property; e.g."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"execution": {
"iopub.execute_input": "2025-05-23T00:19:33.124243Z",
"iopub.status.busy": "2025-05-23T00:19:33.124052Z",
"iopub.status.idle": "2025-05-23T00:19:33.126701Z",
"shell.execute_reply": "2025-05-23T00:19:33.126408Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"TITLE hh.mod squid sodium, potassium, and leak channels\n",
" \n",
"COMMENT\n",
" This is the original Hodgkin-Huxley treatment for the set of sodium, \n",
" potassium, and leakage channels found in the squid giant axon membrane.\n",
" (\"A quantitative description of membrane current and its application \n",
" conduction and excitation in nerve\" J.Physiol. (Lond.) 117:500-544 (1952).)\n",
" Membrane voltage is in absolute mV and has been reversed in polarity\n",
" from the original HH convention and shifted to reflect a resting potential\n",
" of -65 mV.\n",
" Remember to set a squid-appropriate temperature\n",
" (e.g. in HOC: \"celsius=6.3\" or in Python: \"h.celsius=6.3\").\n",
" See squid.hoc for an example of a simulation using this model.\n",
" SW Jaslove 6 March, 1992\n",
"ENDCOMMENT\n",
" \n",
"UNITS {\n",
" (mA) = (milliamp)\n",
" (mV) = (millivolt)\n",
"\t(S) = (siemens)\n",
"}\n",
"\n",
"? interface\n",
"NEURON {\n",
" SUFFIX hh\n",
" REPRESENTS NCIT:C17145 : sodium channel\n",
" REPRESENTS NCIT:C17008 : potassium channel\n",
" USEION na READ ena WRITE ina REPRESENTS CHEBI:29101\n",
" USEION k READ ek WRITE ik REPRESENTS CHEBI:29103\n",
" NONSPECIFIC_CURRENT il\n",
" RANGE gnabar, gkbar, gl, el, gna, gk\n",
" : `GLOBAL minf` will be replaced with `RANGE minf` if CoreNEURON enabled\n",
" RANGE minf, hinf, ninf, mtau, htau, ntau\n",
"\tTHREADSAFE : assigned GLOBALs will be per thread\n",
"}\n",
" \n",
"PARAMETER {\n",
" gnabar = .12 (S/cm2)\t<0,1e9>\n",
" gkbar = .036 (S/cm2)\t<0,1e9>\n",
" gl = .0003 (S/cm2)\t<0,1e9>\n",
" el = -54.3 (mV)\n",
"}\n",
" \n",
"STATE {\n",
" m h n\n",
"}\n",
" \n",
"ASSIGNED {\n",
" v (mV)\n",
" celsius (degC)\n",
" ena (mV)\n",
" ek (mV)\n",
"\n",
"\tgna (S/cm2)\n",
"\tgk (S/cm2)\n",
" ina (mA/cm2)\n",
" ik (mA/cm2)\n",
" il (mA/cm2)\n",
" minf hinf ninf\n",
"\tmtau (ms) htau (ms) ntau (ms)\n",
"}\n",
" \n",
"? currents\n",
"BREAKPOINT {\n",
" SOLVE states METHOD cnexp\n",
" gna = gnabar*m*m*m*h\n",
"\tina = gna*(v - ena)\n",
" gk = gkbar*n*n*n*n\n",
"\tik = gk*(v - ek) \n",
" il = gl*(v - el)\n",
"}\n",
" \n",
" \n",
"INITIAL {\n",
"\trates(v)\n",
"\tm = minf\n",
"\th = hinf\n",
"\tn = ninf\n",
"}\n",
"\n",
"? states\n",
"DERIVATIVE states { \n",
" rates(v)\n",
" m' = (minf-m)/mtau\n",
" h' = (hinf-h)/htau\n",
" n' = (ninf-n)/ntau\n",
"}\n",
" \n",
":LOCAL q10\n",
"\n",
"\n",
"? rates\n",
"PROCEDURE rates(v(mV)) { :Computes rate and other constants at current v.\n",
" :Call once from HOC to initialize inf at resting v.\n",
" LOCAL alpha, beta, sum, q10\n",
" : `TABLE minf` will be replaced with `:TABLE minf` if CoreNEURON enabled)\n",
" :TABLE minf, mtau, hinf, htau, ninf, ntau DEPEND celsius FROM -100 TO 100 WITH 200\n",
"\n",
"UNITSOFF\n",
" q10 = 3^((celsius - 6.3)/10)\n",
" :\"m\" sodium activation system\n",
" alpha = .1 * vtrap(-(v+40),10)\n",
" beta = 4 * exp(-(v+65)/18)\n",
" sum = alpha + beta\n",
"\tmtau = 1/(q10*sum)\n",
" minf = alpha/sum\n",
" :\"h\" sodium inactivation system\n",
" alpha = .07 * exp(-(v+65)/20)\n",
" beta = 1 / (exp(-(v+35)/10) + 1)\n",
" sum = alpha + beta\n",
"\thtau = 1/(q10*sum)\n",
" hinf = alpha/sum\n",
" :\"n\" potassium activation system\n",
" alpha = .01*vtrap(-(v+55),10) \n",
" beta = .125*exp(-(v+65)/80)\n",
"\tsum = alpha + beta\n",
" ntau = 1/(q10*sum)\n",
" ninf = alpha/sum\n",
"}\n",
" \n",
"FUNCTION vtrap(x,y) { :Traps for 0 in denominator of rate eqns.\n",
" if (fabs(x/y) < 1e-6) {\n",
" vtrap = y*(1 - x/y/2)\n",
" }else{\n",
" vtrap = x/(exp(x/y) - 1)\n",
" }\n",
"}\n",
" \n",
"UNITSON\n",
"\n"
]
}
],
"source": [
"print(h.hh.code)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.10"
}
},
"nbformat": 4,
"nbformat_minor": 1
}